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We present the first high angular resolution observations in the near- 
infrared H-b&nd (1.6 /im) of the Luminous Blue Variable star P Cygni. 
We obtained six-telescope interferometric observations with the CHARA 
Array and the MIRC beam combiner. These show that the spatial flux 
distribution is larger than expected for the stellar photosphere. A two 
component model for the star (uniform disk) plus a halo (two-dimensional 
Gaussian) yields an excellent fit of the observations, and we suggest that 
the halo corresponds to flux emitted from the base of the stellar wind. 
This wind component contributes about 45% of the iJ-band flux and 
has an angular FWHM = 0.96 mas, compared to the predicted stellar 
diameter of 0.41 mas. We show several images reconstructed from the 
interferometric visibilities and closure phases, and they indicate a gen- 
erally spherical geometry for the wind. We also obtained near-infrared 
spectrophotometry of P Cygni from which we derive the flux excess com- 
pared to a purely photospheric spectral energy distribution. The if -band 
flux excess matches that from the wind flux fraction derived from the two 
component fits to the interferometry. We find evidence of significant near- 
infrared flux variability over the period from 2006 to 2010 that appears 
similar to the variations in the Ha emission flux from the wind. Fu- 
ture interferometric observations may be capable of recording the spatial 
variations associated with temporal changes in the wind structure. 

Subject headings: stars: variables: other - - stars: early-type - - stars: 
individual (P Cyg, HD 193237) -- stars: winds, outflows -- stars: cir- 
cumstellar matter — stars: mass loss 



1. Introduction 

Luminous Blue Variables (LBVs or S Doradus variables) are evolved, massive 
stars that are characterized by large mas s loss rates and variability on multiple 



timescales. iHumphreys fc Davidson! (11994 ) estimate that the typical mass-loss rate 



of a non-erupting LBV is of the order of M p« 10~ 5 — 1O _3 M yr _1 . The lifetime 
in the LBV phase is uncertain, but is on the order of 25,000 years. One of the 
defining criteria of the LBVs is the detection of a large-scale eruption, when the star 



brightens by several magnitudes. The quiescent times between these eruptions may 
last centuries. In addition to such rare, giant eruptions, these stars also display lesser 
photometric and spec troscopic variations on other timescales, ranging from days to 
decades or centuries (Ivan Genderenll200ll ). The two "prototypical" Galactic LBVs 
are r\ Car and P Cygni, and they probably represent diff erent extremes of both mas s 
and mass loss rate within the scheme of LBV evolution (jlsraelian fc de Grootlll999l ). 



The basi c properties of P Cyg ni (HD 193237, HR 7763, Nova Cyg 1600) were 
estimated by iNajarro et al.l (119971 ) and iNajarrol (120011 ) by comparing ultraviolet, 
optical, and infrared spectroscopic ob servations with results from the non-LTE at- 
mospheric modeling code CMFGEN (IHillier fc Millerl 119981 ). Najarro et al. found 
that P Cyg has a mass- loss rate of 3.2 x 1O~ 5 M yr~ x , a terminal wind speed of 185 
km s _1 , a distance of 1.7 ± 0.1 kpc, and an assumed continuum forming radius of 
75 Rq. The star has a luminosity of (5.6 — 7.0) x 1O 5 L , effective temperature of 
(18.1 — 19.2) kK, and a gravity logoff = 1.20. Their models of the hydrogen and he- 
lium spectral lines in the UV to NIR wavelength range yield a fractional composition 
of nuc/ n R = 0.3, indicat ing that nuclea r processed gas is present in the photosphere. 
The model presented in INajarrol (1200 ll ) also gives estimates of the abundances of C, 
N, O, Mg, Al, Si, Fe, Co, and Ni. Most of these elements have near solar metallicity, 
but C and O are depleted (Nc/N C(S = 0.3 and N /N 0(3 = 0.18) and N is enriched 
relative to solar (N n /N Nq = 6.5). 



Richardson et al.l (120 111 ) examined the optical variabilit y (l^-band photometr y 



and Ha spectroscopy) of P Cyg following the earlier study of iMarkova et al.l (120011 ). 
Both studies found that the Ha emission line increases and decreases in concert with 
the l/-band flux. The variability in the l/-band is due to the portion of the optical 
flux that originates in the stellar win d, but some variab ility may also be due to the 
pulsation properties of the star (e.g., iPercy et al.l 120011 ) . Furthermore, Richardson 
et al. discovered that there were "discrete absorption components" (DACs) in the 
P Cygni absorption trough of the Ha line profile that moved blueward with time, 
indicative of accelerating matter in our line of sight. These were long-lived features 
(2-3 yrs), with a possible recurrence time of ~ 4.7 yrs. This absorption strength 
increases when V is brighter and Ha is stronger. The absorption variations (formed in 
the gas column projected against the star), matched those of the continuum and line 
emission (formed in a large volume surrounding the star) and that was interpreted 
as evidence that the structure of the wind in our direction is comparable to that in 



other directions, i.e., that the wind is spherical in geometry. 



High angular resolution techniques can aid our understanding of the stellar 
winds through measurements of the geometry of the outflows. Stellar winds around 
hot, massive stars are a source of infrared and radio flux from free- free and bound- 
free processes. Because the optical depth of the wind increases with wavelength, 
the observed diameter of the star will appear larger at longer wavelengths (e.g., 
Lamers fc Cassinellilll999l ). The advent of long-baseline optical and near-infrared in- 
terferometry provides a method of directly measuring the extent of emitting regions 
and determining the amount of flux excess at longer wavelengths. Such studies can 
be compared to models of stellar atmos pheres and winds, suc h as those of Hillier & 
Miller (1998; CMFGEN). For example, Iweigelt et all toOl\ ] found evidence of an 
asymmetric wind from the enigmatic LBV rj Car using VLTI/AMBER interferome- 
try. 

Three key studies have used high a n gular resolution techniques to examine the 
Ha emission of P Cygni. IVakili et al.l ( 119971 ) used the Grand Interferometre a 2 
Telescopes to resolve an Ha emitting region of angular diameter 5.5 ± 0.5 mas, which 
corres ponds to 14 i?* for the stellar diameter and di stance deriv e d bvlNaiarro et al. 
(J1997I ). A similar result was recently obtained by iBalan et al.l (120101 ) based upon 
observations made with the Navy Precision Optical Interferometer. They compared 
the angular extent of the Ha emission region to that of the nearby continuum and 
then derived an Ha emitting region that was 3-7 mas in diameter. The Ha emis- 
sio n from the outer region of the wind was explored using adaptive optics techniques 
by IChesneau et al.l (120001 ) , who made a reconstructed image that shows a faint and 
clumpy halo with a radius of ~ 200 mas. The wind struc t ure on the largest scales 
was examined through radio observations by lSkinner et al.l (119971 . Il998l ). Their stud- 
ies resolved circumstellar structures that extend to an angular size of ~ V. The 
structures show a clumpy and/or filamentary appearance. Skinner et al. developed 
a spherically symmetric wind model for the resolved radio nebula, and while they 
were able to obtain a good fit to the spectral energy distribution, their model did 
not explain the structure present in the observations. 

Here we investigate the wind of P Cygni through a comparison of Ha spec- 
troscopy, near infrared H- and iT-band spectrophotometry, and interferometric H- 
band measurements from the CHARA Array using the MIRC beam combiner. We 



outline our spectroscopic measurements in Section 2, and in Section 3 we discuss our 
interferometric measurements and an image reconstruction based upon these data. 
We discuss these results and offer conclusions from our study in Section 4. 



2. Spectroscopy and Photometry 

We collected high resolution Ha spectroscopy and near-infrared spectrophotom- 
etry within 2-3 weeks of our interferometric observations. The Ha spectra were 
obtained from two observatories. The first set of two spectra was obtained at the 
University of Toledo's Ritter O bservatory with the 1 m telescope and echelle spec- 



trograph ([Morrison et al.l 119971 ). The detector was a Spectral Instruments 600 Se- 
ries camera, with a front-illuminated Imager Labs IL-C2004 4100x4096 pixel sensor 
(15 x 15 /iin pixels). Wavelength calibration was accomplished with a Th-Ar discharge 
lamp. These high resolving power (R =26,000) spectra were reduced using standard 
techniques with IRAFi 11 ! with bias frames and flat fields obtained on the same night. 
We reduced three orders of echelle data, which span 6285-6443 A, 6470-6633 A, and 
6666-6834 A. The resulting S/N is roughly 100 in the continuum near Ha. 

A second set of five Ha spectra was obtained at Georgia State University's 
Hard Labor Creek Observatory (HLCO). The data were collected with a 0.5 m RC 
Optical Systems telescope^ and an LHIRES III spectrograph! 13 !. These spectra were 
recorded on a thermo-electrically cooled SBIG-ST8XME CCD. The dispersion was 
accomplished with gratings of either 2400 grooves mm -1 (R ~ 18,000) or 600 grooves 
mm -1 (R ~ 4,500). These data were reduced with standard techniques in IRAF 
utilizing bias, dark, and flat fields. Wavelength calibration was accomplished with a 
built-in Ne discharge lamp in the spectrograph. 

A log of the Ha measurements is listed in Table 1, along with contemporaneous 



11 IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the 
Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the 
National Science Foundation. 



12 http://www. rcopticalsystems.com/telescopes/20truss. html 
13 http://www.shelyak.com/rubrique.php?id_rubrique=6&lang=2 



photoelectric V^-band measurements from the American Association for Variable Star 
Observers (AAVSO). The equivalent widths of the Ha profiles were measured by 
integration between 6 510 and 6617 A for all spectra, which is the same range used by 
Markova et al.l ( 1200 ll ) . Changes in the continuum flux can cause apparent equivalent 
width variations (based upon the changing ratio of emission to continu um flux), so we 



correc ted for the variable continuum in the same manner as done by Markova et al. 



(J200lh and iRichardson et al.l ( 1201 ll ) using the l/-band estimates in Table 1, i.e., 

W x (con) = H/ A (net)10-°- 4(y(t) ~ 4 - 8)) . 



We obtained a single epoch of near-infrared spectrophotometry (Table 2) in the 
K- and L- bands in 2006 with the NASA Inf rared Telescope Facility and SpeX cross- 
dispersed spectrogra ph (jRayner et al.l 120031 ) . This observation (reported earlier by 
Touhami et al.ll201Cf ) was collected with a 3" wide slit. We also obtained five epochs of 
near-infrared spectrophotometry in the H- a nd if -bands using th e Mimir instrument 
on Lowell Observatory's Perkins telescope (jClemens et al.l 120071 ) between 2008 and 
2011. These observations were obtained with a 10" wide slit. These spectra were 
made by combining at least 10 individual spectra obtained by dithering along the 
slit. 

Reductions of the SpeX spectrum were carried out with SpeXtool package 



(jVacca et al.l l2003t ICushing et al. 



2004J ). Reductions of the Mimir data were ac- 



complished using custom software 14 ! that used bias, dark, and flat field frames. We 
corrected for the Mimir detector's non-linearity through a series of flat fields obtained 
on the same observing run to est ablish the response of each pixel to increasing ex- 
posure time (JClemens et al.l 120071 ). The Mimir and SpeX spectra were transformed 
to an ab solute flux scale ( and corrected for telluric absorption) using the xtellcor 
package (JVacca et al.l 120031 ) . This method uses flux calibrator stars of spectral type 
A0 V (in this case HD 192538) that are transformed to flux through reference to a 
model Vega spectrum calculated by R. Kurucz. The transformation takes into ac- 
count rotational broadening, interstellar extinction, and the B and V magnitudes of 
the calibrator (in this case, B = 6.48, V = 6.45). A log of the NIR spectrophotome- 
try is given in Table 2. Note that the spectrophotometry obtained in 2006 and 2008 



14 Available for download at http://people.bu.edu/clemens/mimir/software.html 



was previously published by lTouhami et al.l ( 120101 ). All optical and NIR spectra are 
available upon request. 

The H- and f^-band spectra are plotted in Figure 1, and their rela tive flux place- 
ments confirm the temporal variability found by lTouhami et al.l ( 1201CH ). Figure 1 also 
shows the predicted photospheric SED from a Kurucz model that wa s normalized to 
the b lue part of the spectrum where the wind contributes little flux (jTouhami et al. 
2010). During the course of our observations, the flux excess relative to the photo- 
spheric SED varied between 0.38 and 0.64 mag in the if -band and between 0.63 and 
0.86 mag in the if-band. 

In addition to the optical Ha spectroscopy and the NIR spe ctrophotometry, we 
collect ed \^-band data from the AAVSO and the recent analysis of iPollmann &: Bauer 
((2012J). Figure 2 shows a comparison of the variations over the last six years in 
the l/-band magnitude, Ha emission equivalent width, and IR flux excesses. The 
AAVSO photoele ctric photometry (PEP) m easurements agree well with the \^-band 
measurements of IPollmann fc Bauerl (120121 ). and they show that a local fading oc- 
curred around 2010.0 that was f ollowed by a gradua l incre ase in brightness. Some of 
the y-band measurements from IPollmann fc Bauerl (120121 ) were also reported to the 
AAVSO, so we removed the duplicate points from the AAV SO data set shown in Fig - 
ure 2. The Ha equivalent widt hs from Balan et al. (2010). lRichardson et al.l (120 111 ), 
and IPollmann fc Bauerl (120121 ) are shown together with our new measurements in 
the middle panel (all corrected for the changing continuum). These show many of 
the same trends seen in the \^-band data, and we show a long-term running average 
of the measurements in the top two panels to illustrate this. This was calculated 
using a Gaussian weighting scheme parameterized by a FWHM of 300 days. The 
lower panel shows the H- and K- band flux excesses in the same format, and these 
quantities show evidence that they are correlated with the V^-band magnitudes and 
Ha emission strengths (in particular showing the same fading near 2010.0). These 
trends suggest that all four of these quantities vary in concert accordi n g to c hanges 
in the mass loss rate, as suggested in the analysis of iRichardson et al.l ( 1201 ll ). With 
more NIR data, future analyses may show a direct correlation between the NIR and 
visual magnitudes. 



3. Interferometry 

3.1. Observations 

We obtained interferometric observations of P C yg on three nights, two in 20 10, 
and one in 2011, us ing the MIRC beam combin er (IMonnier et al.l 12004 . 120061 ) at 
the CHARA Array (jten Brummelaar et al.ll2005l ). We observed with MIRC using 
the low-resolution prism (R ~ 42), which disperses the light across eight spectral 
channels in the if-band (1.50—1.75 /im; AA ~ 0.034 /zm for each spectral channel). 
Table 3 presents a log of the observing dates, telescope configurations, range of 
baseline lengths, and observed calibrators. On 2010 Aug 21 and 23, we combined 
the light from three and four telescopes, respectively. On 2011 Sept 3, we combined 
the light from all six telescopes simultaneously. All observations made use of the 
photometric channels installed on MIRC, which measure directly the contribution 
of light from eac h telescope and improve the visibility and closure phase calibration 
Jcheet al.ll2010h . 



To measure the instrument response, we observed calibrator stars with angular 
diameters smaller than 0.7 mas. On each night, we observed a Cyg as the primary 
calibrator; we also used calibrators observed a couple of hours following the P Cyg 
observations to monitor the stability of the visibility calibration. The adopted limb- 
darkened diameters (9r,n) of the calibrators are listed in Table 4. The data were 



reduced using the standard MIRC reduction pipeline (IMonnier et al.ll2007l ). The vis- 
ibilities and closure phases were averaged over the 2—3 min observing blocks. Based 
on an overall assessment of the data quality obtained with MIRC using the photo- 
metric channels, we applied minimum baseline uncertainties of 5% to the squared vis- 



ibilit ies and 0?3 to the closure phases. The calibrated OIFITS data files ( jPauls et al. 
20051 ) are available on request. 



Figure 3 shows the (u, v) coverage on the sky sampled by the CHARA Array 
during the three nights of the P Cyg observations obtained with MIRC. Figure 4 
shows a plot of the squared visibilities measured with MIRC during all three nights. 
The visibilities drop steadily with increasing baseline, and so with spatial frequency, 
indicating that the object is mostly symmetric and resolved on the longest baselines. 
Figure 5 shows the closure phases measured on each closed triangle. There are small 
non-zero closure phases (~ 2°) on some triangles, possibly indicating the presence of 



a small asymmetry in the light distribution. 



3.2. Geometric Models 

Our results in Figure 2 show that the star's Ha emission was stronger during 
the 2011 CHARA observations than during the 2010 CHARA observation. However, 
the visibilities measured by CHARA/MIRC for both sets in the if -band are com- 
parable. Therefore, in addition to fitting the data from each epoch separately, we 
also combined the two CHARA data sets to constrain better the models. As Figure 
6 shows, we fit three types of geometrical models to the interferometric visiblities: 
a single uniform disk, a single circularly symmetric Gaussian, and a two-component 
model where the star is represented as a uniform disk and the extended wind emission 
is represented by a circular Gaussian centered on the star. In the two-component 
model, we fi xed the stellar disk d iame ter at flup = .411 mas (75 R Q at 1.7 kpc), 



as found by iNajarro et al.l (119971 ) and iNajarrol (1200 ll ) , and solved for the FWHM 



size of the Gaussian wind (^fwhm) and the flux ratio between the wind and the star 

(, /wind / /star ) ■ 

Table 5 lists the results for the 2010 and 2011 epochs separately as well as for 
the combined data set. In Figure 4, the comparison of the models to the visibilities 
demonstrates the superiority of the two-component model. In both epochs as well as 
the combined dataset, the xl is significantly improved in the two-component model 
compared with a single uniform disk or a single circular Gaussian. We attempted 
to fit an elliptical Gaussian to the data as well, but we found that the \ 2 of the 
fit was not improved significantly. Moreover, the difference between the FWHM of 
the major and minor axes was less than 5%, with a deviation of less than 1.2 cr, 
suggesting that the wind is essentially circular in the if -band continuum. 

Figure 6 shows images of the best fit models for the uniform disk, circular 
Gaussian, and the two-component model for the combined data set. The best fit two- 
component model for the combined 2010 and 2011 data sets gives a wind FWHM size 
of 0.96 ± 0.02 mas [R = 2.3i?*), contributing approximately 45% of the total flux. 
The parameter uncertainties quoted in Table 5 include the small spread in the results 
due to the uncertainties in the assumed photospheric angular diameter (A^tjd/^ud ~ 
3%). Our geometric model assumes that the circular Gaussian representing the wind 
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is centered on the star. The star may in fact block flux from the wind behind the star. 
To account for this effect, we calculated the total flux in the wind that is coincident 
with the star (/ cen t)- We corrected the flux contributions by removing 0.5/ cent from 
the wind component and adding 0.5/ ccnt to the stellar flux (to roughly account for 
the foreground/background portions of the wind). After applying this correction, we 
find that the wind contribution drops to 42% of the total flux in the combined data 
set. 

We attempted to fit for the uniform disk diameter of the star as a free parameter 
in the two-component model, but the results were not consistent between epochs, 
with values ranging from 0.39 to 0.62 mas. The resolution of the CHARA Array on 
the longest baseline in the if -band is ~ 0.52 mas, therefore, the stellar diameter is 
only marginally resolved and would require a more precise calibration of the inter- 
ferometric visibilities to measure reliably. The situation is further complicated by 
model degeneracies between the stellar diameter, size of the wind, and flux contri- 
bution of each component. For instance, assuming a larger stellar angular diameter 
results in a model that contains a fainter, but more extended wind. Because we fixed 
the stellar diameter over a narrower range of parameter space than is allowed for by 
the interferometeric data alone, we suspect that the uncertainties listed in Table 5 
do not fully represent the model degeneracies. Therefore, differences in the param- 
eters between the two epochs should not be treated as significant. Additionally, the 
sampling of the (u, v) coverage may affect the parameters of the fit. We note that 
the reduced xl of the combined fit is close to 1, indicating that the model does a 
good job reproducing the data from each epoch. 



3.3. Image Reconstructions 

Figure 5 shows the closure phases measured with the MIRC beam combiner at 
the CHARA Array. The small, but non-zero, closure phases show that there may 
be some structure in the wind. The geometric models computed in Section 3.2 are 
point-symmetric and do not account for the non-zero closure phases. We attempted 
to map the asymm etry by reconstruct ing images of P Cygni using the MArkov Chain 



IMager fM ACIM; llreland et al. 



20061 ) and BiSpectrum Maximum Entropy Method 



(BSMEM; iBaron fc Youngl 120081 ) software packages. MACIM randomly moves flux 
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within a pixel grid to reconstruct the image; the movement is regulated by a simulated 
annealing temperature. BSMEM uses a gradient descent method to converge to the 
best image. 

Examples of different image reconstructions of P Cygni based on the MIRC 
data from 2010—2011 are shown in Figure 7. We combined the epochs to maximize 
the (u, v) coverage, but also reconstructed images for each epoch separately to avoid 
potentially blurring the motion of structures within the wind. The top row of Figure 7 
shows a comparison of images reconstructed using BSMEM and MACIM. For the 
BSMEM reconstruction, we started with the initial image and a prior set to a 2.0 
mas Gaussian. Note that the initial image defines the starting position of the flux 
distribution, while the prior image defines the probability of where the flux is likely 
to move during the reconstruction. For the MACIM reconstruction, we used our best 
fit uniform disk and Gaussian model as the initial image. In both of these cases, we 
find a larger amount of flux in the central region of the image and a more compact 
size for the extended emission, compared to the two-component geometric models. 
This could suggest that the boundary between the star and the wind is more blurred 
than in our simpler models. However, this could also be the result of the software 
having difficulty reconstructing a sharper boundary between the "edge" of the star 
and the fainter, more diffuse emission. The faint, extended tails in the north-south 
and east-west directions line up roughly with gaps in the (u, v) coverage (compare 
the images with combined coverage for 2010—2011 in Fig. 3) and are most likely 
artifacts produced in the reconstruction process. The generation of such artifacts 
could also be influenced by small, baseline-dependent calibration errors, such as the 
visibilities near B/X ~ (50 — 70) x 10 6 in Figure 4, which are systematically below 
the model fit. 

The MACIM software allows for simultaneously fitting for a uniform disk while 
reconstructing the extended emission. Including a uniform disk with a diameter of 
0.411 mas that contributes 55% of the flux, MACIM produced the images in the 
lower panels of Figure 7. For both reconstructions, we used the Gaussian component 
of our two-component model as the initial image for the extended emission. MACIM 
offers the flexibility of using different regularizers that reduce the weighting for non- 
physical images by balancing the tasks of lowering the \ 2 statistic while optimizing 
the regularization statistic. The image on the lower left of Figure 7 used the MACIM 
option of a "compressed sensing regularizer," which minimizes spatial gradients in 
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the image (e.g., iDonohol 120061 ; ICandes fc Tad 120061 ) . For the image on the lower 
right, we used our Gaussian model as a prior to define the probability for where the 
emission might be located, in order to keep the emission more centrally located. In 
both cases, the sharp edge between the star and the wind is retained because we 
included a fixed uniform disk in the reconstruction process. These images show that 
the if -band emission from the wind of P Cygni is largely spherical and consistent 
with the overall size and shape that we derived from our geometric model. 

In Figure 5, we overplotted the closure phases computed from the image in the 
lower-left panel of Figure 7. All of the image reconstructions, for the combined data 
set as well as the individual epochs, provide reasonable fits to the small, non-zero 
closure phases, but unfortunately, we could not find a unique solution to describe the 
location and shape of a slight asymmetry in the wind. The nearly point-symmetric 
morphology of P Cygni (with closure phases close to 0°) makes it difficult to constrain 
the detailed structure of the wind. Additionally, diffuse emission from an extended, 
Gaussian envelope is a difficult case for interferometric imaging, because the image 
reconstruction techniques can generate artifacts associated with gaps in the frequency 
coverage. 



Vakili et al.l (119971 ) examined P Cygni interferometrically using the Grand In- 



terferometre a 2 Telescopes (GI2T) in the Ha and He I A 6678 lines. They deduced 
that there was a structure in the wind located at a projected radial separation of 
R fa 4R* (0.8 mas) away from the star. The angular resolution (FWHM) of the 
CHARA Array with the MIRC beam combiner is « 0.5 mas, so we should detect 
such a structure if it is relatively bright compared to the star an d wind. The clo- 
sure phase (in radians) is ~ -^asymmetric /-^symmetric ( IMonnierl 120071 ). and our results 
(see Fig. 5) show that the largest closure phase we measure is ~ 2°, or 0.035 radi- 
ans. If we assume the entire closure phase quantity is due to a single asymmetry in 
the wind (such as a blob or clump), then any such blob would contribute less than 
about 4% of the symmetric stellar and wind flux. According to our models of the 
wind halo (Fig. 8), the wind flux at R = 0.8 mas is quite faint in the if -band (only 
about 1% of the maximum light). Consequently, observations like ours would only 
detect rather extreme and isolated wind density enhancements (i.e., a distribution 
of clumps surrounding the star would yield a smaller net closure phase because the 
flux distribution would appear more symmetrical). 
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3.4. Limits on the Presence of a Binary Companion 



Kashil (j2010l ) argued that 17th-Century eruptions of P Cygni might have resulted 
from an interaction with a B-type, main-sequence st ar in a 7- year, hi ghly elliptical or- 
bit. However, the historical light curve reported by lde Grootl ( 119881 ) . which was used 
as a basis in these models, has been re-evalu ated and has a more "typical" appear- 
ance of an LBV eruption (jSmith et al.ll201ll ) when viewed with a sparsely sampled 
light curve. While no periodic radial velocity variation was found in the recent high - 
resolution spectroscopic analysis over a 15-year period by iRichardson et al.l (120 ll[ ). 
the possibility of a binary is not elimi nated given the high incidence of massive stars 
in binary systems ( jMason et al.ll2009i ). A companion in a 7 yr orbit would probably 
have an angular semimajor axis of approximately 6 mas for the probable distance 
and mass of P Cygni, and a binary with such a separation might be detectable in 
our interferometric observations. 

To examine that possibility, the high-precision closure phases and visibilities 
measured with the MIRC beam combiner at the CHARA Array were evaluated to 
place limits on the presence of a binary companion to P Cygni. A binary star will 
produce a periodic signature in the visibilities and the closure phases, where the 
frequency of the variation depends o n the separation of the com ponents and the 
amplitude depends on their flux ratio (JBodenll2000l ; iMonnierl 120071 ) . We focused our 
efforts on the data set from 2011 Sept 3, to avoid the motion of the hypothetical 
companion between 2010 and 2011. Additionally, the 2011 data offers better (u,v) 
coverage and more closure phase triangles for computing the detection limits. 

We investigated two possible scenarios. In both cases, we fixed the diameter 
of the primary stellar component to be 0.41 mas (75 R Q ) and assumed that the 
secondary is a point source. In the first scenario, we explored whether the small, 
non-zero closure phases could be accounted for by a binary companion alone. We 
used our two-component uniform disk and symmetric Gaussian model of the wind 
emission optimized for the 2011 epoch (6*fwhm = 0.898 mas). We then fit for a 
binary system by searching through a grid of separations in RA and Dec. over a 
range of ±14 mas and solving for the secondary-to-primary flux ratio of the binary 
at each step. We also allowed the flux contributions from the primary star and its 
wind to vary in order to accommodate the additional flux from the hypothetical 
secondary. We ran through two iterations of the grid search: on the first pass we 
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used a step size of 0.2 mas in separation and allowed the RA and Dec. separations 
to vary to their best fit values; on the second pass we used a fixed step size of 0.01 
mas to finely map the x 2 surface near the location of the absolute minimum. Using 
this approach, we found a best-fit binary solution where the agreement with the 
visibilities is similar to the symmetric two-component model and the y 2 calculated 
from the closure phases is reduced from 198.4 for the symmetric model, where we 
do not fit for the closure phases, to Xcp = 99.5 when including a binary companion 
(based on 160 closure phase measurements). However, the image reconstructions for 
this epoch give a x 2 calculated from the closure phases of only ~ 54.8. Therefore, 
while a binary companion could account for some of the non-zero closure phase 
signal, the image reconstructions that map the fine-scale structure in the wind do a 
better job of fitting the data. Based on the analysis of the binary fits with a two- 
dimensional Gaussian wind, we estimate that any possible binary companion must 
be more than 4.9 mag fainter than the central star in P Cygni or 5.6 mag fainter 
than the star+wind combined. 

In the second scenario, we selected a MACIM image reconstruction of the wind 
assuming a uniform disk central star for the data from 2011 Sept 3 and investigated 
whether adding a binary component would improve the fit. We added in a binary 
model to the image of the wind and searched through a grid of separations in RA 
and Dec. over a range of ±14 mas and solved for the secondary-to-primary flux ratio 
of the binary at each step while allowing the flux contribution from the wind to vary. 
We found a best-fit binary solution where the total x 2 calculated from the visibilities 
and closure phases is reduced from 154.3 (x y 2 = 99.5, Xcp = 54.8) for the image+UD 
to x 2 — 145.5 (xy2 = 93.6, Xcp = 51.9) when including a binary companion (based 
on 120 visibilities and 160 closure phase measurements). Performing an F-test on 
the ratio of the reduced xl values (0.527/0.553 = 0.95), we find that this ratio can be 
exceeded by about 35% of random observations, suggesting that the improvement in 
the fit by adding the binary parameters is not significant. Based on the analysis of 
the binary fits using the reconstructed image of the wind, we find a tighter restriction 
on the presence of a binary companion in that it must be more than 5.3 mag fainter 
than the central star or 6.0 mag fainter than the star+wind combined in the if-band. 

We can estimate the absolute magnitude of the LBV star in P Cygni and then de- 
termine limits on the kinds of faint companions that remain undetected. Based upon 



the if-band flux measurements given in Table 2 and the calibration of I Cohen et al. 
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(120031 ). we estimate that the apparent i?-band magnitude was 3.28 ± 0.09 during 
2010 to 2011. Then, given the distance and reddening from iNajarro et al.l (119971 ). 
the absolute H-band magnitude at that time was —8.14 ± 0.15 for the LBV and 
its wind, or —7.49 ± 0.15 for the LBV alone (based upon the flux fraction from the 
results for the two-component model in Table 5). The AH = 5.3 mag limit from 
above then suggests that we would have probably detected any main sequence star 
brighter than if = —2.2 mag. This magnitude corresponds approximately to that of 
a Bl V star ( iCoxl 120001 ) . Thus, our results appear to rule out main-sequence com- 
panions of types O to B0V, but not later. In all cases, a secondary star very close 
(within a few i?*) to the primary would blend with the primary and would be unde- 
tectable. Future studies with more CHARA/MIRC data utilizing all six telescopes 
will produce stronger constraints on the abse nce of a co mpanion. The number of 
known LBVs wit h a ste ll ar com panion is small (jVinkll2012l ). While this analysis does 
not dispr ove the iKashil ( 20101) conjecture, the improved light curve of the eruption 
shown bvlSmith et al.l (20111 ). the lack of periodic radial velocity variation found by 
Richardson et al.l ( 1201 ll ). and this analysis point toward a single-star nature for P 
Cygni. 



Discussion and Conclusions 



Our CHARA Array observations provide the first high angular resolution look at 
the LBV P Cygni in the if -band continuum. We find that the angular size of the wind 
is much larger than the 0.41 mas diam eter predicted for the photosphere fr om spectral 
models and the distance of the star ( Najarro et al.lll997l ; lNajarroll200ll ). A spatial 
flux model consisting of a uniform disk photosphere and a circular Gaussian halo 
provides an excellent match to the interferometric observations. The halo probably 
corresponds to flux emitted in the base of the stellar wind of P Cygni. The FWHM 
of the halo light is approximately 1 mas, which is smaller than the 5.5 mas size found 
for the Ha emission (form ed over a larger wind volume because of its higher optical 
depth; IVakili et al.lll997h . 



The two-component model also provides an estimate of the relative flux contri- 
butions of the star (uniform disk) and wind (circularly symmetric Gaussian), and 
the flux ratio of 1.36:1 (all epochs) from interferometry is consistent with that de- 
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rived from spectrophotometry. The interferometry obtained in 2010 yields a ra- 
tio of F stSLT /F wind = 1.64 ± 0.09. This implies a predicted magnitude difference 
of AH = 2.51og[(F wind + i r s tar)/-^star] = 0.52 ± 0.05 mag. Our spectrophotome- 
try obtained within a month of the 2010 interferometric data shows an excess of 
AH = 0.55 ±0.01 mag (Table 2), in excellent agreement with the interferometric re- 
sults. Thus, the SED models used to determine the IR flux excess and the geometric 
models used in the interferometric analysis agree within the uncertainties. 

Our interferometric results obtained with the MIRC instrument and the CHARA 
Array include precise closure phases (Fig. 5). These measurements reveal the pos- 
sibility of a small amount of asymmetry present in the wind and may help probe 
the geometry of the stellar wind outflow. From image reconstructions, we found 
there may be some wind asymmetry (as indicated by small non-zero closure phases), 
but the wind is spherical within observational limits. It may be possible in future 
studies, with more complete coverage of the (u, v) plane than we obtained (Fig. 3), 
to determine the exact location of any bright asymmetries in the wind. The image 
reconstruction process could also be aided by the development of a regularizer that 
would keep the image smooth in azimuth in order to avoid the problem of diffuse, 
extended emission from mapping onto gaps in the (u, v ) coverage. 

Mo dels of LBV atmospheres and winds can a ccurately reproduce the emergent 



spectra (INaiarro et al.lll997l : iHillier & Millerl ll998l. We utilized the CMFGEN model 



of P Cygni with the parameters of lNajarrol ( 200ll ) to compare the theoretical visibility 



curve with the observed visibility curve (Fig. 4). We computed the visibility curve for 
the derived wind radial light distribution in the if-band continuum for the adopted 
mass loss rate M = 3.2 x 10~ 5 M o yr~ 1 , a stellar radius of 76 R Q and a distance of 1.7 
kpc. The CMFGEN model predicts a radial profile with a half width half maximum 
of R/R & = 88 (Fig. 8). However, the visibilities associated with the model made a 
poor match with the observed values, so we rescaled the angular size until it best 
matched the visibilities (Fig. 4). This resulted in the CMFGEN model being scaled 
17% larger than expected. We overplot the visibilities of the scaled CMFGEN model 
as the dotted red line in Figure 4, and show the scaled model for comparison in 
Figure 6. The resulting agreement has a reduced xt °f 1-5? only marginally worse 
than that of our two-component model (xt = 1-1) ■ The major deviation in the scaled 
curve is at the largest baselines, where we sample the smallest angular scales. We 
also used this scaled CMFGEN model as an initial image for image reconstruction, 
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and we found that the resulting image was nearly identical to that which used our 
two-component model as the initial image. 

We compare the calculated radial light distribution from the CMFGEN model 
to our radial distribution from the image reconstructions in Figure 8. The radial 
profile from our reconstructed images includes a sharp edge from the star. In reality, 
the sharp edge would be smoother due to effects of electron scattering and free-free 
emission. Therefore, in order to build a more realistic profile, we performed a Gaus- 
sian smoothing to the images. The smoothed profiles appear similar to the predicted 
profile from the CMFGEN models but the scales of the inner light distribution differ. 
However, we caution that the angular resolution of our interferometry is insufficient 
to differentiate between the models at scales < 0.5 mas. The general agreement at 
large scales shows that the CMFGEN model prediction about the spatial distribution 
of the wind flux is consistent with our observations. 

The model radial intensity distribution had to be rescaled to fit the visibilities, 
which resulted in a size about 17% larger than the predicted size. There are several 
plausible explanations for this difference. First, P Cygni could be closer to us than the 
assumed distance of 1.7 kpc, so that its angular size appears larger. Secondly, there 
are a large number of hydrogen emission lines in the if -band (see Fig. 1) that form in 
the wind at larger radii than the continuum flux. The emission lines are blended with 
the continuum flux in our low spectral resolution observations with MIRC (R ~ 42), 
but their net contribution may lead to an overestimate of the size of the spatial 
intensity distribution compared to what would be observed for the continuum alone. 
Finally, it is possible that some adjustments to the wind model could account for the 
difference in the angular size of the wind flux (if the a bove reasons are insuffic ient). 



The mass-loss rate of P Cygni is known to vary (e.g., iRichardson et al.ll201ll ). and 
the observed Ha strength was high during the second epoch of CHARA observations 
(Fig. 2). This might indicate that the actual mass loss rate was larger than the 
assumed model value. In addition, small changes related to the assumed velocity 
law and/or wind clumping factor could also lead to differences in the angular size 
prediction comparable to the observed difference. 

Our results represent the first images of the circumstellar environs close to the 
prototypical LBV, P Cygni. The wind appears to be spherically symmetric. The 
results of the interferometric analysis and the spectroscopic analysis are mutually 



consistent, so long term temporal variations in wind emission should also be detected 
in spatial observations through interferometry. P Cygni should remain a prime tar- 
get for monitoring with optical and NIR spectroscopy, photometry, and improved 
interferometry. 
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Table 1. Ha Spectroscopy 







Date 




w x 




w x 




UT 




(HJD-2,400,000) 




(net) 


V 


(corr) 




(YYYY-MM-DD) 


Source 


(d) 


R 


(A) 


(mag) 


(A) 




(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 




2008-10-22 


Ritter 


54762 


26,000 


-87.4 


4.72 


-94.1 


to 
to 


2010-08-29 


Ritter 


55438 


26,000 


-85.8 


4.81 


-85.0 




2011-05-10 


HLCO 


55692 


4,500 


-89.6 


4.67 


-101.0 




2011-05-24 


HLCO 


55706 


18,000 


-86.7 


4.67 


-97.7 




2011-06-14 


HLCO 


55727 


4,500 


-86.2 


4.62 


-101.7 




2011-06-21 


HLCO 


55734 


4,500 


-84.1 


4.58 


-103.0 




2011-08-16 


HLCO 


55790 


4,500 


-89.6 


4.66 


-101.9 





Table 2. NIR Spectrophotometry 



UT 




Date 


log(F H ) 


AH 


log(F K ) 




AK 


Date 


Observatory/ 


(HJD-2,400,000) 


(1.629 urn) 


(1.629 Mm) 


(2.179 fim) 


(2 


.179 urn) 


(YYYY-MM-DD) 


Instrument 


(d) 


(ergs s _1 cm -2 A -1 ) 


(mag) 


(ergs s" 1 cm -2 A -1 ) 




(mag) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 




(7) 


2006-09-16 


IRTF/SpcX 


53994 






-11.62 




0.72 g 


2008-10-20 


Lowcll/Mimir 


54759 


-11.21 


0.64 


-11.57 




0.86 


2009-07-13 


Lowcll/Mimir 


55025 


-11.27 


0.48 


-11.62 




0.73 


2009-11-06 


Lowell/Mimir 


55141 


-11.31 


0.38 


-11.66 




0.63 


2010-07-03 


Lowcll/Mimir 


55380 


-11.25 


0.55 


-11.59 




0.81 


2010-11-27 


Lowell/Mimir 


55527 


-11.25 


0.55 


-11.60 




0.80 




MODEL 




-11.46 


0.00 


-11.91 




0.00 



Table 3. CHARA/MIRC Observations 



UT Date 


Telescope 


Min. 


Baseline 


Max. 


Baseline 


Calibrator 




(YYYY-MM-DD) 


Designations 




(m) 




(m) 


Names 


to 


2010-08-20 


S1-W1-W2 




108 




279 


a Cyg, C Cas 




2010-08-23 


S2-E2-W1-W2 




108 




251 


a Cyg, 7 And, ( Cas 




2011-09-03 


S1-S2-E1-E2-W1-W2 




34 




331 


a Cyg, 7 And, 9 Cas 
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Table 4. Adopted Angular Diameters for the Calibrators 



Name HD 



(mas) 



Reference 



o Cyg 


202850 


0.574 ± 0.017 


1 




7 And 


219080 


0.659 ± 0.017 


2. 


3,4,5 


( Cas 


3360 


0.307 ± 0.021 


6 




flCas 


6961 


0.608 ± 0.019 


7 





Referen ces. - - (1) ISchaefer et all J20iqh : (2) 



Che et al. (20 



Bonneau et al- 



ii): (3) iBarnes et all (I1978h: (4) 



2006h: ( 5 ) iKervella fc Fouaue 



( 2008 ): (6) IWeigelt et all ( 120071 ): (7) iGies et al. 
fJ2007t ). 
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Table 5. Model Fit Results 



Parameter 


UT 2010 Aug 20+23 UT 2011 Sep 3 


Combined Fit 




Uniform Disk Model 




6» UD (mas) 


0.902 ± 0.007 0.838 ± 0.004 


0.858 ± 0.005 


/Cu 


6.21 9.28 


9.04 




Circularly Symmetric Gaussian Model 




#fwhm (mas) 


0.558 ± 0.004 0.524 ± 0.003 


0.536 ± 0.003 


A/y 


4.57 5.83 


5.99 


Two-Component Model: Uniform Disk and Circular Gaussian 


/wind 


0.396 ± 0.014 0.475 ± 0.020 


0.450 ± 0.018 


f a 
Jcorr, wind 


0.379 0.443 


0.423 


#fwhm (mas) 


1.121 ± 0.027 0.898 ± 0.022 


0.964 ± 0.022 


/star 


0.604 ± 0.014 0.525 ± 0.020 


0.550 ± 0.018 


f a 

,/coit, star 


0.621 0.557 


0.577 


6» UD (mas) 


0.41 (fixed) 0.41 (fixed) 


0.41 (fixed) 


A// 


0.69 1.13 


1.13 


Num. Vis. 


120 120 


240 



a These values correspond to the corrected flux where flux from a circular 
Gaussian wind behind the star would be blocked from our line of sight, so half 
of the flux in the Gaussian coincident with the star is subtracted from the wind 
and added to the star. 
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Fig. 1. — if- (left) and K- band (right) spectrophotometry. Dates are given in the 
legend and match the observations listed in Table 2. The (black) spectrum of low 
flux is the Kurucz model used to determine the IR excess. A color version is available 
in the online edition. 
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Fig. 2. — \^-band photometry from the AAVSO and iPollmann fc Bauer! (120121 ) is 
shown in the top panel (note that many of the measurements of Pollmann & Bauer 
were also reported to the AAVSO, making the data sets have some overlap, so only 
those measurements from Pollmann & Bauer are shown in such cases), Ha equivalent 
width (corrected for a changing continuum) in the middle panel, and the HK-ba.nd 
IR excesses in the bottom panel. We over plotted a running average of the V^-band 
flux (red dotted line) in both the top and middle panel, as well as a similar curve (blue 
dashed line) for the Ha measurements in the second panel to show the similarity of 
the variability of these measurements. The epochs of CHARA/MIRC observations 
are marked with vertical lines in the bottom panel. A color version is available in 
the online edition. 
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Fig. 3. — (u,v) coverage during the MIRC observations of P Cyg for the 2010 and 
2011 data, as well as for the combined data set. 
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Fig. 4. — Calibrated visibilities measured for P Cyg using MIRC on UT 2010 Aug 
20+23 and 2011 Sep 3. The black dashed-dotted line is the best fit uniform disk 
model, the blue dashed line represents the best Gaussian model, the solid green line 
is our two-component model, and the dotted red line is the rescaled CMFGEN model 
(see discussion). A full color version is available in the online version. 
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Fig. 5. — Closure phases measured for P Cyg using MIRC on UT 2010 Aug 20+23 
and 2011 Sep 3. The solid lines show the closure phases associated with the recon- 
structed image in the lower left panel of Figure 7. On 2010 Aug 23, we obtained 
two observations on P Cyg separated in time by about 30 minutes. There is a small 
amount of rotation in the (u, v) plane between these observations which samples the 
reconstructed image in a slightly different way. The two solid lines in the plots for 
this date show the reconstructed closure phases for each data point. On 2010 Aug 
22 and 2011 Sep 3 we only obtained one observation of P Cyg. 
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Fig. 6. — Model flux distributions of P Cyg fit to the interfero metric visibilities in the 
combined data from 2010-2011 : a uniform disk (upper left), a circularly symmetric 
Gaussian (upper right), and a two-component model (bottom left). In the bottom 
right we show the flux distribution of the CMFGEN model scaled to optimize the 
fit to the interferometric visibilities (the visibility curve shown in Fig. 4). Contour 
intervals are drawn at 0.01%, 0.05%, 0.1%, 0.5%, 1%, 2%, 4%, 6%, 8%, 10%, 20%, 
30%, 40%, and 50% of the peak flux in each panel (contours extend up to 60% and 
70% of the peak flux for the circular Gaussian and 60%, 70%, and 80% of the peak 
flux for the CMFGEN model). 
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Fig. 7. — Image reconstructions of P Cygni based on our MIRC data from 2010 and 
2011 (240 visibility measurements and 232 closure phase measurements). Top left: 
Image reconstruction from BSMEM using a 2.0 mas Gaussian as the initial image 
and prior (xv 2 = 331.4 for the visibilities, Xcp = ^0-4 f° r the closure phases, where 
the x 2 is calculated as the the number of observations minus the degree of the fit). 
Top right: Image reconstruction from MACIM using our two-component geometric 
model as the initial image (xy2 = 267.1, Xcp = 82.1). Bottom left: MACIM image 
reconstruction of the extended emission while fitting for a uniform disk of 0.41 mas 
that contributes 55% of the light (x y 2 = 218.7, Xcp = 86.1). This reconstruction 
used the Gaussian component of our two component model as the initial image and 
used a regularizer to minimize spatial gradients in the reconstructed image. Bottom 
right: MACIM image reconstruction of the extended emission made assuming a 
stellar flux component with a UD of 0.41 mas that contributes 55% of the light (xy2 
= 194.9, Xcp = 103.1). This reconstruction used the Gaussian component from our 
two-component model as the initial image and as a prior to define the probability for 
how the flux moves during the reconstruction process. Contour intervals are drawn 
at 0.05%, 0.1%, 0.5%, 1%, 2%, 3%, 4%, 5%, 6%, 7%, 8%, 9%, 10%, 20%, 40%, 60% 
of the peak flux in each panel. In all cases, the x 2 calculated from the closure phases 
is smaller than the number of measurements. However, we suspect that this is the 
result of the closure phases being so close to and that small movements in the flux 
during the reconstruction process can reproduce the signal in many different ways, 
allowing the software to find a very precise, but not necessarily reliable, solution. 
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Fig. 8. — Theoretical light distribution of the star and wind in the if -band from 
CMFGEN models (red solid curve) directly from the CMFGEN model, as well as 
a rescaled version from a fit to the visibility curve (red dotted curve). The aver- 
age radial profile of the MACIM image reconstruction (lower left panel of Fig. 7) 
is shown by the diamond symbols, along with Gaussian-smoothed versions of the 
reconstruction over 10, 20, and 30% of the stellar radius (blue dashed lines). The 
purple dot-dashed curve represents the two-component model derived from the vis- 
ibilities. Finally, the radial distribution of the flux from the image reconstruction 
with no regularizer and no uniform disk component (top right image of Fig. 7) is 
plotted as plus signs. 



